Dispersion interactions from a local polarizability model 
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A local approximation for dynamic polarizability leads to a nonlocal functional for the long- 
range dispersion interaction energy via an imaginary-frequency integral. We analyze several local 
polarizability approximations and argue that the form underlying the construction of our recent van 
der Waals functional [O. A. Vydrov and T. Van Voorhis, Phys. Rev. Lett. 103, 063004 (2009)] 
is particularly well physically justified. Using this improved formula, we compute dynamic dipole 
polarizabilities and van der Waals Ce coefficients for a set of atoms and molecules. Good agreement 
with the benchmark values is obtained in most cases. 



I. INTRODUCTION 

Recently we developed [l| a nonlocal correlation en- 
ergy functional that describes the entire range of van der 
Waals interactions in a general and seamless fashion, us- 
ing only the electron density and its gradient as input. 
Improving upon its predecessors 0, Q , the new van der 
Waals density functional denoted VV09, has a simple 
analytic form, generalized to spin-polarized systems and 
well-behaved in some important limits. In the asymp- 
totic long-range regime, VV09 reduces to a form similar 
to the models of Refs. [1] and yet with some crucial 
differences. In this article, we examine this long-range be- 
havior in detail and present some test results of dynamic 
dipole polarizabilities and asymptotic van der Waals Cg 
coefficients. 



II. FORMALISM 

For two compact systems A and B separated by a large 
distance R, the nonretarded dispersion interaction energy 
@ behaves asymptotically as —Cq^R~^ with the Cg co- 
efficient given by the formula 



(1) 



where a{iu) is the average (isotropic) dynamic dipole po- 
larizability at imaginary frequency iu. A simple but often 
sufficiently accurate approximation is to describe a by a 
local model: 



(2) 



a{iu) = / dra(r,m). 



The long-range dispersion interaction energy between 
systems A and B can then be written Q in terms 
of local polarizabilities as 
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1,1, , a(r,iu) a(r' ,iu) , . 
dii / dr / dr'— — '—^ — (3) 



where r is within the domain of system A and r' is within 
the domain of B. 



In Refs. 0, 11 i, a simple model for a(r, iu) was 
derived from the response properties of a uniform elec- 
tron gas (UEG). The zero wave vector UEG dielectric 
function at frequency w is given by 



1 _ ^ 



(4) 



where lo^ — Aiine^ / m is the plasma frequency for the 
electron density n. In nonuniform systems, the local ana- 
log of ujp can be defined via Wp(r) — Ann{Y)e'^ Jm. Then 
the local polarizability for uj = iu is. found as 0, [13] 
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(5) 



Plugging Eq. ([5|) into Eq. ([3]) we arrive at the Andersson- 
Langrcth-Lundqvist (ALL) formula [3| 
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An immediately apparent problem with Eq. ([5]) is its 
treatment of static polarizability: 
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dr a(r, 0) = / dr 
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a{v, 0) is constant everywhere, therefore the above inte- 
gral is divergent unless a cutoff is introduced. Eq. 
taken as it is, yields finite but severely overestimated 
-Edisp- These difficulties are circumvented d, [l^, [HI 
by the introduction of sharp density-based integration 
cutoffs in Eqs. ^ and Calculated polarizabilities 

and Cq coefficients are admittedly 0, [3 

sensitive to the 

choice of the cutoff criterion, although the prescription of 
Refs. H. [l0l.[Til appears to work well in many cases. Note 
that Ref. [Ifll gave separate cutoff criteria for the spin- 
compensated and the fully spin-polarized cases. To our 
knowledge, a prescription for a general spin-polarization 
case has never been put forth. 

An integration cutoff discards density tail regions, 
which is not entirely satisfactory from the formal point 
of view. In the u — >■ (X) limit, the /-sum rule requires 
that 



-r ^ 

a{iu) -> 
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where N is the number of electrons in the system. Omis- 
sion of the density tails leads to the reduction of N, in 
violation of the /-sum rule. 

Another formal shortcoming of the ALL theory was 
pointed out by Nesbet [T3l.[l3|. who argued that a more 
appropriate relation between a and e is given by the 
Clausius-Mossotti formula 



1 
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so that Eq. ^ should be replaced by 



a{r,iu) = — 



47r a;^(r)/3 + u"^ 



(9) 
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It appears that Nesbet 's articles went unnoticed, because 
in the numerous practical applications [l5l - [3l| of the ALL 
formula, Nesbet's suggestion was never utilized. 

The validity of Eq. (flUl) is corroborated by the example 
of interacting jellium spheres. For two identical spheres 
of uniform density and radius tq separated by distance 
R (such that i? 3> tq) the interaction energy is given by 



(11) 



The above result is exactly reproduced if Nesbet's model 
of Eq. ((TU]) is plugged into Eq. whereas the ALL 
formula ^ underestimates this result by the factor of 
3-\/3 « 5. We note in passing that all three versions of 
the vdW-DF functional of Refs. fail to reproduce 

Eq. PT|) even on the qualitative level, yielding incorrect 
dependence on the electron density. 

Local polarizability models of Eqs. (O and ([TU)) were 
derived using the UEG dielectric function of Eq. 
UEG is rather dissimilar to our target systems — 
molecules. UEG has a continuous excitation spectrum 
and a zero band gap (i.e. it is a metal) , whereas molecules 
have a discrete spectrum with a gap between the ground 
state and the fist excited state. The polarizability model 
could be made more realistic by introducing a gap. For a 
semiconductor with a band gap hujg , the zero wave vector 
dielectric function [34] is typically written as 



e(a;) = 1 + 



9 7 ' 



(12) 



Using this e{uj) in the Clausius-Mossotti formula 1^, we 
obtain 
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c-p'(r) 



(13) 



where we introduced a "local gap" hujg{r). The above 
a(r, iu) leads via Eq. ([3]) to the energy expression 
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where loq = 



aj2/3. A suitably chosen ujg{r) obvi- 



ates any need for an integration cutoff in E'disp and a(iu). 
As a result, the /-sum rule on a{iu) is obeyed. 

An apt model for ujg{r) can be deduced by examining 
the behavior of the electron density n(r). In atoms, n(r) 
can be approximated as piecewise^exponential. In the 
density tails, the exact behavior [35] is known: 



n(r) ~ exp (— Q;|r|) , with 



(15) 



where / is the ionization potential. Generalizing the re- 
sult of E iq. (ITS I. we can define a "local ionization poten- 
tial" as |3£ 
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Taking Uujg{r) oc /(r), in Ref. pLj] we defined 
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where C is an adjustable parameter. We fitted C to a 
benchmark set of 17 van der Waals Cg coefficients and ob- 
tained f]\ the optimal value of C = 0.0089. It is instruc- 
tive to consider the ratio Wg(r)//(r) = 8V0.0089 = 0.755. 
This ratio seems reasonable since ujg should be somewhat 
smaller than /. In the uniform density limit, Eq. (jl7p 
gives ojg = 0, so that our a(r,m) of Eq. (fT3|) reduces to 
Nesbet's a{r,iu) of Eq. (fTO]) . 

Eqs. (HH), (Ull), and dTT]) require only the total elec- 
tron density as input and include no dependence on spin- 
polarization. The question of the proper treatment of 
spin, not fully resolved in the ALL theory [l3|, does not 
arise in this model. We mention in passing that vdW- 
DF functionals of Refs. 0, d, were defined only for 
the spin-compensated case and their extension to spin- 
polarized systems is nontrivial. In VV09 the depen- 
dence on spin-polarization enters only at shorter range. 
In the long-range limit, the VV09 nonlocal correlation 
energy reduces to Eq. (|14p . In this regard, a clarifica- 
tion should be made: The coefficient before the double 
integral in Eq. (|H)) is twice the coefficient in Eq. (13) of 
Ref. [ij because these formulas compute different things. 
Eq. (|14p computes the interaction energy between sys- 
tems A and B, hence the integral over r is limited to the 
part of space confining system A, while the integral over 
r' is limited to the domain of B. On the other hand, 
Eq. (13) of Ref. [1] gives the nonlocal correlation energy, 
which includes inter- and intramolecular contributions, 
hence both r and r' integrals are over the entire space. 



III. BENCHMARK TESTS 

All calculations reported in this section were performed 
at the LC-wPBE08 0| electron densities (using oj = 
(14) 0.45 Uq^, as suggested in Ref. 42]), except for the H 
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TABLE I. Isotropic polarizabilities a{iu) calculated via 
Eqs . ((2| and (I13[) compared to the reference values for BeH2 
[4(j| and BH [43|. Atomic units are used. 



BeHa BH 



u 


Ref. 


Calc. 


Ref. 


Calc. 


0.0 


19.760 


19.413 


21.430 


17.442 


0.142857 


17.234 


17.049 


16.582 


15.506 


0.333333 


11.280 


11.509 


10.343 


10.738 


0.6 


6.084 


6.440 


5.625 


6.139 


1.0 


2.917 


3.161 


2.755 


3.044 


1.666667 


1.233 


1.360 


1.198 


1.307 


3.0 


0.418 


0.480 


0.424 


0.453 


7.0 


0.081 


0.105 


0.091 


0.096 



atom polarizability, computed at the Hartree-Fock (i.e. 
exact in this case) density. For the numerical integration, 
we use the Euler-Maclaurin-Lebedev unpruned (75,302) 
quadrature grid. The aug-cc-pVQZ basis set is used in 
all calculations. All the numbers in this section are given 
in atomic units (a.u.). 

Using the VV09 model, given by Eqs. ^ and (IT71) . 
we have calculated the isotropic dynamic dipole polariz- 



abilities as functions of imaginary frequencies for several 
atoms and small molecules for which accurate reference 
data t38..-.41.j are available. For LiH, BeH2, and BH, we 
used the same bond lengths as in Refs. [H, li^,!!^. The 
results are given in Fig.[T]and TablelH The agreement be- 
tween the calculated and reference values of a{iu) is gen- 
erally quite good, although this method has a tendency of 
underestimating static polarizabilities a(0). The largest 
errors in a(0) are observed for LiH (Fig. [T|) and for 
alkali- metal atoms (not shown). Underestimation of a(0) 
causes rather large errors in Cg coefficients for alkali- 
metal atoms, as shown below. 

When the distance between species A and B is large 
compared to the size of these systems, |r — r'|~^ in 
Eq. can be taken out of the integral as R~^, leading 
to the -C^^R^^ form with C^^ given by Eq. To 
further assess the quality of the VV09 local polarizability 
model, we have calculated isotropic dispersion Cg coeffi- 
cients for a number of atoms and molecules. As expected 
from Eq. (U]), any errors in the polarizability a"^{iu) are 
reflected in Cq^ and similarly in Cq^. It is sufficient 
to include only in our benchmark set, since the ac- 
curacy for C^'^ and Cg^ determines the accuracy for 
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TABLE II. Cq^ coefficients (a.u.) for closed-shell species calculated by several methods. Experimental geometries [43| are used 
for all molecules. MPE stands for the mean percentage error and MAPE stands for the mean absolute percentage error. 



Molecule 


vdW-DF-04^ 


vdW-DF-09'' 


vdW-DF-10'= 


vvog"^ 


Accurate 


Ref.'= 


He 


2.93 


1.63 


0.76 


1.45 


1.46 


[39] 


Ne 


9.45 


6.52 


3.07 


8.44 


6.35 


m 


Ar 


62.67 


61.41 


25.29 


70.08 


64.42 


m 


Kr 


114.3 


120.0 


47.7 


131.2 


130.1 


[44] 


Be 


269 


330 


102 


186 


214 


[39] 


Mg 


649 


835 


246 


425 


627 


[39] 


Zn 


269 


240 


87 


163 


284 


[45] 


Ha 


16.82 


12.53 


5.09 


10.28 


12.09 


[46] 


Na 


78.76 


77.59 


31.96 


88.70 


73.43 


[46] 


Cla 


289.3 


336.8 


131.4 


366.7 


389.2 


[47] 


HP 


23.12 


18.01 


7.97 


21.13 


19.00 


[48] 


HCl 


114.3 


119.9 


47.2 


124.6 


130.4 


[48] 


HBr 


180.1 


198.2 


76.1 


200.2 


216.6 


[48] 


CO 


87.56 


86.34 


35.01 


93.51 


81.40 


[49] 


COa 


127.6 


130.6 


54.5 


159.4 


158.7 


[49] 


CSa 


586.3 


731.7 


274.3 


739.4 


871.1 


[50] 


ocs 


316.8 


370.1 


143.4 


395.6 


402.2 


[50] 


NaO 


136.1 


140.3 


58.4 


172.4 


184.9 


[51J 


CH4 


122.0 


130.1 


50.8 


129.6 


129.6 


[52] 


CCI4 


1436 


1882 


715 


2044 


2024 


[53] 


NH3 


82.47 


79.32 


32.00 


82.78 


89.03 




HaO 


46.96 


40.83 


17.17 


44.95 


45.29 


M 


SiH4 


338.0 


406.1 


147.2 


344.6 


343.9 


M 


SiP4 


360.9 


382.7 


158.6 


455.8 


330.2 


[55] 


HaS 


186.6 


208.9 


79.1 


200.3 


216.8 


[56] 


SO2 


239.5 


265.1 


106.5 


305.2 


294.0 


[50] 


SFe 


568.0 


659.7 


274.6 


869.9 


585.8 


[57] 


C2H3 


191.3 


210.3 


81.0 


210.3 


204.1 


[58] 


U2rl4 




zyo.o 




007 Q 

zy / .0 


Qfin 
oUU.z 


159{ 


CaHe 


330.4 


386.1 


148.8 


396.6 


381.8 


[52] 


CH3OH 


194.0 


208.5 


83.2 


226.1 


222.0 


[60] 


CH3OCH3 


458.7 


532.5 


207.7 


567.9 


534.0 


m 


Cyclopropane 


480.7 


596.1 


228.3 


632.6 


630.8 


[52] 


CeHg 


1297 


1715 


647 


1838 


1723 


[58] 


MPE (%) 


-2.8 


-0.5 


-60.9 


1.2 






MAPE (%) 


18.5 


10.4 


60.9 


10.7 







^ The method of Ref. J. 
The method of Ref. @|. 

This method is denoted as vdW-DF2 in Ref. [3I . 
^ The formalism proposed in Ref. [J| and described in this work. 
Literature references for the accurate benchmark Cg^'* values. 



Cg For example, VV09 strongly underestimates the 
Cg coefficient for the Li-Li interaction, and as a result, 
all Cg coefficients for Li interacting with other species 
are also underestimated. On the contrary, VV09 gives 
very accurate Cg coefficients for He-He and Kr-Kr, and 
consequently, Cg for He-Kr is also very accurate. 

In Ref. [1] we reported the C^^ coefficients for a set 
of 17 closed-shell species, computed within the VV09 
methodology. In fact, the value of C = 0.0089 in Eq. ([17]) 
was fitted to that set. In this study, we test whether this 
fit is transferable to atoms and molecules outside of the 
training set. In Table ITll we assembled a set of 34 closed- 
shell species for which accurate C^^ are known. Using 
this benchmark set, we compare the accuracy of VV09 



[T| to the similar methods of Refs. 0, S H^- In the 
asymptotic limit, all these methods reduce to the form of 
Eq. ()14l) . but with different models for wq, as discussed 
in Ref. 3] . Deviations from the reference values are sum- 
marized in Table [III as mean (signed) percentage errors 
(MPE) and mean absolute percentage errors (MAPE). 
VV09 and vdW-DF-09 exhibit very similar accuracy with 
MAPE of just over 10%. vdW-DF-04 is somewhat less 
accurate with MAPE of 18.5%. The latest reparame- 
terization (denoted as vdW-DF2 in Ref. [Ill, but called 
vdW-DF-10 here for consistency) yields very poor Cg co- 
efficients: as compared to the reference values, vdW-DF- 
10 underestimates C^^ by a factor of 2.6 on average. 

The good performance of vdW-DF-09 for Cg coeffi- 
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TABLE III. Ce ^ coefficients (a.u.) for open-shell species cal- 
culated using Eq. (|14p . 



Molecule 


VV09 


Accurate 


Ref. 


H 


6.75 


6.50 


[39] 


Li 


565 


1389 


[39] 


B 


87.6 


99.5 


[45] 


C 


47.0 


46.6 


[45] 


N 


27.65 


24.10 


[51J 


O 


18.19 


14.89 


[51] 


F 


12.21 


9.52 


[45] 


Na 


669 


1556 


[39] 


Al 


353 


528 


[45] 


Si 


253 


305 


[45] 


P 


179 


185 


[45] 


S 


130 


134 


[45] 


CI 


94.7 


94.6 


[45] 


Ga 


255 


498 


[45] 


(je 


zoi 


Q C /I 

354 


|45J 


As 


222 


246 


[45] 


Se 


190 


210 


[45] 


Br 


158 


162 


[45] 


O2 


66.18 


61.57 


[5^ 


NO 


77.83 


69.73 


[5^ 


MPE (%) 
MAPE (%) 


-9.8 
18.7 







cients motivated Sato and Nakai [62] to devise a pairwise 
atom-atom dispersion correction using the local polariz- 
ability model 3] underlying the construction of vdW-DF- 
09. We believe that the VV09 model of Eq. ^ can also 
be successfully employed in this scheme. 

As mentioned above, none of the three versions of 
vdW-DF 0, 0, Issj has been generalized for open-shell 
systems, whereas VV09 is defined for a general spin- 
polarized case. In Table IIIIl we compare the C^^ co- 



efficients predicted by VV09 to the accurate reference 
values for 20 open-shell species. The agreement is sat- 
isfactory in most cases. The largest errors are observed 
for the alkali-metal atoms Li and Na. The strong under- 
estimation of the Cg coefficients for alkali-metal atoms 
was also noted for the ALL formula It is likely that 
the local approximation of Eq. ^ is inadequate for such 
highly polarizable systems as alkali metals. 

IV. CONCLUSIONS 

The ALL formula © for the long-range dispersion 
energy enjoys growing popularity il^-i3|| , even though 
it has been superseded by more general [l|-y, [SJl and 
more accurate (63l . [63 | methods. A simple change from 
Eq. ([6]) to Eq. p4)) improves the theory in several impor- 
tant ways: the sharp integration cutoff is obviated and 
consequently the /-sum rule is recovered; the model sys- 
tem of two distant jellium spheres is properly described; 
accurate Cg coefficients are predicted for many atoms 
and molecules including open-shell species. Eq. de- 
scribes the asymptotic limit and has to be damped at 
short range. To this end, empirical damping functions 
are often used (see e.g. Ref. ^2]). 

The general and seamless van der Waals functional 
VV09 lf| reduces to Eg. (|14p in the large separation limit. 
As our recent study [63 shows, VV09 performs well not 
only in the asymptotic limit, but also near equilibrium 
intermonomer separations, provided that an adequate ex- 
change functional is used. 
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